Guión - Prueba Didáctica COA 15087-05

Contexto

La presente clase forma parte de un concurso de oposición abierto de la UNAM, en el cual se evaluará la asignatura optativa Propiedades Físicas de los Sólidos (clave: 0087).

Dicha asignatura tienen como objetivo introducir a los estudiantes a los aspectos formales de las propiedades formales de los sólidos desde el punto de vista de la teoría cuántica haciendo uso de los conceptos desarrollados en el curso obligatorio de Química Cuántica I.

En esta clase no puedo presuponer que ustedes ya llevaron dicho curso, por lo que intentaré hacer mención de los resultados conocidos en lugar de realizar deducciones, posiblemente laboriosas.

La duración aproximada de la clase es de 20 minutos.

El tema a exponer y evaluar es el Unidad 5/5.3 Método de bandas, o más generalmente conocido como la Teoría de Bandas. Dicha teoría explica que los estados electrónicos en estructuras cristalinas generan bandas y su estructura nos permite describir el comportamiento de los electrones. Es un hito en el estudio de los sólidos que introdujo Felix Bloch en 1929 al aplicar la teoría cuántica a los sólidos [ Über die Quantenmechanik der Elektronen in Kristallgittern].

Objetivo de la clase

Describir cómo surge una banda electrónica con lo cual podremos definir las bandas electrónicas, y decir qué estructura tienen. Ocuparemos la estructura de las bandas para explorar las propiedades ópticas.

1 Planteamiento

El punto de partida en todo problema de la mecánica cuántica es escribir el Hamiltoniano (el Hamiltoniano es un operador que describe la energética del sistema), si el Hamiltoniano no depende explícitamente el tiempo las soluciones de la ecuación de Schrödinger

\[ \begin{align} i\hbar \frac{\partial }{\partial t} \Psi = \hat{H}\Psi \end{align} \]

pueden separarse en la variable temporal y la variable espacial de la siguiente forma

\[ \begin{align} \Psi(\vec{r},t) = e^{iEt/\hbar} \psi(\vec{r}) \end{align} \]

donde \(\psi\) es solución de la ecuación de Schrödinger estacionaria

\[ \begin{align} \hat{H} \psi = E \psi \end{align} \]

Por lo que, el problema a resolver es, dado un \(\hat{H}\), determinar quienes son las funciones que cumplen que al aplicarle el operador \(\hat{H}\) me regresan la misma función multiplicada por un factor. Las funciones se llaman eigenfunciones, \(\psi\), y los factores que multiplican, \(E\) se conocen como eigenvalores.

Comencemos escribiendo el Hamiltoniano de un cristal. Un cristal se conforma por electrones, neutrones y protones. En nuestro contexto, estamos interesados en describir la fenomenología de los movimientos colectivos más que la descripción individual de los entes que constituyen al cristal; por lo que, consideraremos los electrones de valencia y a cada núcleo con los electrones de las capas llenas como una partícula, al cual llamaremos core. Entonces, el Hamiltoniano se escribe como

\[ \begin{align} \hat{H} &= \sum_{\text{cores}} (\text{Energía cinética})_{\text{core}}+\sum_{\text{electrones}} (\text{Energía cinética})_{\text{el}} \\ &+\sum_{\text{core,core}} (\text{Energía potencial})_{\text{core,core}} \\ &+ \sum_{\text{core,el}} (\text{Energía potencial})_{\text{core,el}}\\ &+ \sum_{\text{el,el}} (\text{Energía potencial})_{\text{el,el}} \end{align} \]

Recordemos que la energía potencial entre dos cargas puntuales \(q_1\) y \(q_2\) separados una distancia \(r_{12}\) es de la forma,

\[ \begin{align} U_{12} = \frac{1}{4\pi\epsilon_0} \frac{q_1q_2}{r_{12}} \end{align} \]

Técnicamente, si denotamos las posiciones de los core por \(\vec{R}_i\) y su momento por \(\vec{P}_i\), y la posición de los electrones por \(\vec{r}_i\) y su momento por \(\vec{p}_i\) el Hamiltoniano de \(N\) cores y \(n\) electrones se escribe de la siguiente forma

\[ \begin{align} \hat{H} &= \sum_{i=1}^{N} \frac{\hat{P}_i^2}{2m_{\text{core},i}} + \sum_{i=1}^{n} \frac{\hat{p}_i^2}{2m_{\text{e}}}+\sum_{i=j+1}^{N}\sum_{j=1}^{N-1} \frac{1}{4\pi\epsilon_0} \frac{q_{\text{core},i}q_{\text{core},j}}{|\vec{R}_i-\vec{R}_j|}\\ &+\sum_{i=1}^{N}\sum_{j=1}^{n} \frac{1}{4\pi\epsilon_0} \frac{q_{\text{core},i}(-e)}{|\vec{R}_i-\vec{r}_j|}+\sum_{i=j+1}^{N}\sum_{j=1}^{N-1} \frac{1}{4\pi\epsilon_0} \frac{(-e)(-e)}{|\vec{r}_i-\vec{r}_j|} \end{align} \]

Nuestro primer paso es realizar la aproximación de Bohr-Oppenheimer, la cual permite la separación del movimiento de los electrones del movimiento de los núcleos. Es decir, debido a las diferencias de masa, el tiempo característico de los electrones es mucho menor que el tiempo característico de los núcleos. Por lo que los núcleos están estáticos para los electrones y los núcleos ven a los electrones como la contribución promedio de los electrones.

En dicha aproximación, el Hamiltoniano asociado a la mecánica de los electrones es:

\[ \begin{align} \hat{H}\approx \hat{H}_{\rm BO} &= \sum_{\text{electrones}} (\text{Energía cinética})_{\text{el}} \\ &+ \sum_{\text{core,el}} (\text{Energía potencial})_{\text{core,el}}\\ &+ \sum_{\text{el,el}} (\text{Energía potencial})_{\text{el,el}} \\ &+\sum_{\text{core,core}} (\text{Energía potencial})_{\text{core,core}} \end{align} \]

donde el último término de \((\text{Energía potencial})_{\text{core,core}}\) es constante.

El usual reescribir al Hamiltoniano de Bohr-Oppenheimer de la siguiente manera:

\[ \begin{align} \hat{H}_{\rm BO} &= \hat{H}_{\text{el}} +\sum_{\text{core,core}} (\text{Energía potencial})_{\text{core,core}}\\ &= \hat{H}_{\text{el}} + U_{\text{core,core}} \end{align} \]

Como se planteo al inicio, nuestro objetivo es resolver,

\[ \begin{align} \hat{H}_{\text{BO}} \psi = E \psi\\ (\hat{H}_{\text{el}} + U_{\text{core,core}}) \psi = E \psi\\ \hat{H}_{\text{el}} \psi + U_{\text{core,core}} \psi = E \psi \\ \hat{H}_{\text{el}} \psi = (E-U_{\text{core,core}}) \psi \\ \hat{H}_{\text{el}} \psi = E_{\text{el}} \psi \\ \end{align} \]

donde hemos definido \(E_{\text{el}}=E-U_{\text{core,core}}\).

En consecuencia, nuestro objetivo ahora es resolver \(\hat{H}_{\text{el}} \psi = E_{\text{el}} \psi\) , el cuál explícitamente es

\[ \begin{align} \hat{H}_{\text{el}} &= \sum_{i=1}^{n} \frac{\hat{p}_i^2}{2m_{\text{e}}}\\ &+\sum_{i=1}^{N}\sum_{j=1}^{n} \frac{1}{4\pi\epsilon_0} \frac{q_{\text{core},i}(-e)}{|\vec{R}_i-\vec{r}_j|}\\ &+\sum_{i=j+1}^{N}\sum_{j=1}^{N-1} \frac{1}{4\pi\epsilon_0} \frac{(-e)(-e)}{|\vec{r}_i-\vec{r}_j|} \end{align} \]

El último término es difícil de tratar, por lo que despreciaremos en primera aproximación dicho término. Entonces,

\[ \begin{align} \hat{H}_{\text{el}} &\approx \sum_{i=1}^{n} \frac{\hat{p}_i^2}{2m_{\text{e}}}\\ &+\sum_{i=1}^{n}\sum_{j=1}^{N} \frac{1}{4\pi\epsilon_0} \frac{q_{\text{core},j}(-e)}{|\vec{R}_j-\vec{r}_i|} \end{align} \]

Notamos que podemos reescribir el Hamiltoniano como suma individual de Hamiltonianos de una sola partícula.

\[ \begin{align} \hat{H}_{\text{el}} &\approx \sum_{i=1}^n \bigg[ \frac{\hat{p}_i^2}{2m_{\text{e}}} - \sum_{j=1}^{N} \frac{1}{4\pi\epsilon_0} \frac{eq_{\text{core},j}}{|\vec{R}_j-\vec{r}_i|} \bigg]\\ \hat{H}_{\text{el}} &\approx \sum_{i=1}^n \bigg[ \hat{H}_{\text{el}, i} \bigg] \end{align} \]

Otro resultado que utilizaremos de la teoría cuántica es que, las soluciones del Hamiltoniano \(\hat{H}_{\text{el}}\) que se puede escribir como la suma de Hamiltonianos de partículas no interactuantes es

\[ \begin{align} \psi(\vec{r}_1,\ldots,\vec{r}_n) &= \displaystyle \prod_{i=1}^n \psi_i(\vec{r}_i)\quad \text{y}\\ E_{\text{el}} &= \displaystyle \sum_{i=1}^n E_{\text{el},i} \end{align} \]

Notamos que el \(\hat{H}_{\text{el},i}\) es igual para todos.

Idependientemente del formalismo, el resultado es que si resolvemos el Hamiltoniano de un electrón \(\hat{H}_{\text{el},i}\) entonces tenemos la solución del Hamiltoniano de \(n\) electrones \(\hat{H}_{\text{el}}\).

Por lo que ahora nuestro problema es resolver

\[ \hat{H}_{\text{el},i} = \frac{\hat{p}_i^2}{2m_{\text{e}}} - \sum_{j=1}^{N} \frac{1}{4\pi\epsilon_0} \frac{eq_{\text{core},j}}{|\vec{R}_j-\vec{r}_i|} \]

Ahora ya no es necesario hacer la distinción con el subíndice \(i\) es decir, podemo reescribir

\[ \begin{align} \hat{p}_i \to \hat{p}\\ \hat{r}_i \to \vec{r} \end{align} \]

Entonces, \[ \begin{align} \hat{H}_{\text{el},i} = \frac{\hat{p}^2}{2m_{\text{e}}} - \sum_{j=1}^{N} \frac{1}{4\pi\epsilon_0} \frac{eq_{\text{core},j}}{|\vec{R}_j-\vec{r}|} \end{align} \]

Identificamos que el Hamiltoniano en el esquema de la función de onda es, \[ \begin{align} \hat{H}_{\text{el},i} = -\frac{\hbar^2}{2m_{\text{e}}}\nabla^2 - \sum_{j=1}^{N} \frac{1}{4\pi\epsilon_0} \frac{eq_{\text{core},j}}{|\vec{R}_j-\vec{r}|} \end{align} \tag{1}\]

del cual buscamos las soluciones,

\[ \begin{align} \hat{H}_{\text{el},i} \psi(\vec{r}) = E \psi(\vec{r}) \end{align} \tag{2}\]

1.1 Ejemplo: Molecula \({\rm H}_2\)

En el caso de la molécula de \({\rm H}_2\) el Hamiltoniano a resolver,

\[ \begin{align} \hat{H}_{\text{H}_2} = -\frac{\hbar^2}{2m_{\text{e}}}\nabla^2 - \frac{1}{4\pi\epsilon_0} \frac{e^2}{|\vec{R}_1-\vec{r}|} - \frac{1}{4\pi\epsilon_0} \frac{e^2}{|\vec{R}_2-\vec{r}|} \end{align} \tag{3}\]

Notamos que si tuvieramos sólo los primeros dos términos del Hamiltoniano, tenemos un Hamiltoniano de un átomo de Hidrógeno en la posición \(\vec{R}_1\), y la solución para el nivel de energía más bajo sería el \(1s\). Entonces,

\[ \bigg( -\frac{\hbar^2}{2m_{\text{e}}}\nabla^2 - \frac{1}{4\pi\epsilon_0} \frac{e^2}{|\vec{R}_1-\vec{r}|}\bigg) \phi_{1s}(\vec{r}-\vec{R}_1) = E_{1s} \phi_{1s}(\vec{r}-\vec{R}_1) \]

De igual forma, tomado el primer término de Equation 3 y el último, la solución de más baja energía es un orbital \(1s\) ahora en la posición \(\vec{R}_2\)

\[ \bigg( -\frac{\hbar^2}{2m_{\text{e}}}\nabla^2 - \frac{1}{4\pi\epsilon_0} \frac{e^2}{|\vec{R}_2-\vec{r}|}\bigg) \phi_{1s}(\vec{r}-\vec{R}_2) = E_{1s} \phi_{1s}(\vec{r}-\vec{R}_2) \]

La gráfica de ambos se muestra en Figure 1.

Code

from pylab import *
from scipy.special   import sph_harm
from scipy.integrate import nquad
from scipy.special   import factorial
from scipy.special   import hyp1f1
import plotly.graph_objects as go


π = pi
ħ = 1.054571817e-34 # J.s = kg m^2 /s^2 * s = kg m^2 /s
m = 9.10938356e-31  # kg
ħ2_2m = ħ*ħ/2/m    # kg^2 m^4 /s^2 /kg = kg m^4 /s^2
# print (ħ2_2m,"J m^2")
ħ2_2m = ħ2_2m*(6.242e18)*(1e20)
# print (ħ2_2m,"eV Å^2")
e = 1.602e-19          # C
ɛ0 = 8.8541878176e-12  # C^2/(N m^2)
e2_4πɛ0 = e*e/4/π/ɛ0
# print (e2_4πɛ0, "J m")
e2_4πɛ0 = -e2_4πɛ0*(6.242e18)*(1e10)
# print (e2_4πɛ0, "eV Å")

def realYlm(θ,ϕ,l=0,m=0):
    # https://scipython.com/blog/visualizing-the-real-forms-of-the-spherical-harmonics/
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html
    Y = sph_harm(abs(m),l,ϕ,θ)
    if m<0:
        Y = sqrt(2) * (-1)**m * Y.imag
    elif m>0:
        Y = sqrt(2) * (-1)**m * Y.real
    return Y

def R(r,n=1,l=0,Znuclear=1,a0=5.29177e-1):
    Zn = 2*Znuclear/n/a0
    ρ  = 2*Znuclear*r/(n*a0)
    N  = 1/factorial(2*l+1)*sqrt( Zn*Zn*Zn*factorial(n+l)/(2*n*factorial(n-l-1)) )
    R  = ρ**l*exp(-ρ/2)*N*hyp1f1(l+1-n,2*l+2,ρ)
    # https://en.wikipedia.org/wiki/Confluent_hypergeometric_function
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.hyp1f1.html
    return R

def Ψnlmxyz(n,l,m,X,Y,Z,X0=0,Y0=0,Z0=0):
    r = sqrt((X-X0)**2+(Y-Y0)**2+(Z-Z0)**2)
    θ = arccos((Z-Z0)/r)
    ϕ = arctan2((Y-Y0),(X-X0))+π
    return R(r,n=n,l=l)*realYlm(θ,ϕ,l=l,m=m)

N = 100*1J
xmin,xmax = -10,10
ymin,ymax = -10,10
zmin,zmax = -10,10

X,Y = mgrid[xmin:xmax:N,ymin:ymax:N]
Z = zeros_like(X)

# Geometría
H2 = { 0:array([-0.50000, 0.00000,0.00000]),
       1:array([+0.50000, 0.00000,0.00000])}
# """

PLOTS = []
for key in H2:
    ni,li,mi = 1,0,0
    sitioRi  = H2[key]
    Xi,Yi,Zi = sitioRi
    Ψnlm_i   = Ψnlmxyz(ni,li,mi,X,Y,Z,X0=Xi,Y0=Yi,Z0=Zi)
    PLOTS.append(go.Surface(x=X,y=Y,z=Ψnlm_i.real,opacity=0.5,showscale=False))
    PLOTS.append(go.Scatter3d(x=[Xi],y=[Yi],z=[0],mode='markers',showlegend=False,marker=dict(color='red')))

fig = go.Figure( data=PLOTS )

camera = dict(
    eye=dict(x=0.5, y=2, z=0.1)
)

fig.update_layout(
    width=500, height=500,
    scene = dict(
        xaxis_title = 'x',
        yaxis_title = 'y',
        zaxis_title = '1s',
        xaxis = dict(range=[-2,2],showticklabels=False),
        yaxis = dict(range=[-2,2],showticklabels=False),
        zaxis = dict(nticks=2,showticklabels=False)
    ),
    scene_camera=camera
)
fig.show()

Figure 1: Gráfica de los orbitales 1s centrados en \(\vec{R}_1\) y \(\vec{R}_2\).

Por tanto, podemos proponer la solución de \(\hat{H}_{\text{H}_2}\)

\[ \begin{align} \psi =c_1 \phi_1 + c_2 \phi_2 \end{align} \]

1.2 Ejemplo: Cuadrado de Hidrógenos (modelo de juguete)

Code

from pylab import *
matplotlib.rcParams['text.usetex'] = True

# Geometría
H2 = { 0:array([-0.50000, -0.50000, 0.00000]),
       1:array([+0.50000, -0.50000, 0.00000]),
       2:array([+0.50000, +0.50000, 0.00000]),
       3:array([-0.50000, +0.50000, 0.00000])
       }
# """

fig,ax = plt.subplots(figsize=(4,4))
eps = 0.08
for key in H2:
    ax.scatter( H2[key][0],H2[key][1],c='#cd7687',s=200,zorder=2 )
    ax.text(  H2[key][0]+eps,H2[key][1]-eps, 
    r'${\rm H}_{(}$'+'$_{0}$'.format(key)+r'$_)$',fontsize=22 )

for (i,j) in [ (0,1),(1,2),(2,3),(3,0) ]:
    ax.plot( [H2[i][0],H2[j][0]],[H2[i][1],H2[j][1]],c='black',lw=6,zorder=1,alpha=0.8 )
ax.set_axis_off()

\[ \begin{align} \mathcal{H} = \begin{pmatrix} E_0 & -t & 0 & -t \\ -t & E_0 & -t & 0 \\ 0 & -t & E_0 & -t \\ -t & 0 &-t & E_0 \\ \end{pmatrix} \end{align} \]

Code

from pylab import *

E0 = -27
t = 20
N = 4
D = ones(N)*(E0)
L = ones(N-1)*(-t)
H = diag(D) + diag(L,1)+diag(L,-1)
H[(0,N-1)] = -t
H[(N-1,0)] = -t
E = eigvalsh(H)

fig,ax = plt.subplots(figsize=(4,4))
for n in range(len(E)):
  ax.plot([0,1],[E[n],E[n]],c='red',alpha=0.5)
ax.text(-0.4,E0,"$E_0$",fontsize=20)
ax.text(-0.4,E0+2*t,"$E_0+2t$",fontsize=20)
ax.text(-0.4,E0-2*t,"$E_0-2t$",fontsize=20)
ax.set_axis_off()
ax.set_xlim(-1,2);

1.3 Ejemplo: Cuadrado de Hidrógenos (modelo de juguete)

Code

from pylab import *

E0 = -27
t = 20
N = 100
D = ones(N)*(E0)
L = ones(N-1)*(-t)
H = diag(D) + diag(L,1)+diag(L,-1)
H[(0,N-1)] = -t
H[(N-1,0)] = -t
E = eigvalsh(H)

fig,ax = plt.subplots(figsize=(4,4))
for n in range(len(E)):
  ax.plot([0,1],[E[n],E[n]],c='red',alpha=0.5)
ax.text(-0.4,E0,"$E_0$",fontsize=20)
ax.text(-0.4,E0+2*t,"$E_0+2t$",fontsize=20)
ax.text(-0.4,E0-2*t,"$E_0-2t$",fontsize=20)
ax.set_axis_off()
ax.set_xlim(-1,2);

1.4 Método de variación lineal (\(\text{H}_2\))

Buscamos resolver,

\[ \begin{align} \hat{H}_{\text{H}_2} \psi = E \psi \end{align} \]

En lugar de resolverlos analíticamente el problema, podemos utilizar el método varicional lineal. Proponemos una solución que sea combinación lineal de \(N\) funciones conocidas

\[ \begin{align} \psi = \sum_{i=1}^{N} c_i \phi_i \end{align} \]

y buscamos qué coeficientes minimizan la energía. Esto nos lleva a un sistema de ecuaciones homegéneas

\[ \begin{align} \mathcal{H} \mathcal{C} = E\mathcal{S} \mathcal{C} \end{align} \]

donde la matriz

\[ \begin{split} \mathcal{H} = \begin{pmatrix} H_{11} & H_{12}\\ H_{21} & H_{22} \end{pmatrix} \end{split} \]

\[ \begin{split} \mathcal{S} = \begin{pmatrix} S_{11} & S_{12}\\ S_{21} & S_{22} \end{pmatrix} \end{split} \]

Proponemos que la solución sea, \[ \begin{align} \psi = c_1\phi_1 + c_2\phi_2 \end{align} \]

Donde \(\phi_1(\vec{r}-\vec{R}_1)\) y \(\phi_2(\vec{r}-\vec{R}_2)\) son funciones \(1s\) localizadas en \(\vec{R}_1\) y \(\vec{R}_2\).

Code

from pylab import *
from scipy.special   import sph_harm
from scipy.integrate import nquad
from scipy.special   import factorial
from scipy.special   import hyp1f1


π = pi
ħ = 1.054571817e-34 # J.s = kg m^2 /s^2 * s = kg m^2 /s
m = 9.10938356e-31  # kg
ħ2_2m = ħ*ħ/2/m    # kg^2 m^4 /s^2 /kg = kg m^4 /s^2
# print (ħ2_2m,"J m^2")
ħ2_2m = ħ2_2m*(6.242e18)*(1e20)
# print (ħ2_2m,"eV Å^2")
e = 1.602e-19          # C
ɛ0 = 8.8541878176e-12  # C^2/(N m^2)
e2_4πɛ0 = e*e/4/π/ɛ0
# print (e2_4πɛ0, "J m")
e2_4πɛ0 = -e2_4πɛ0*(6.242e18)*(1e10)
# print (e2_4πɛ0, "eV Å")

def realYlm(θ,ϕ,l=0,m=0):
    # https://scipython.com/blog/visualizing-the-real-forms-of-the-spherical-harmonics/
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html
    Y = sph_harm(abs(m),l,ϕ,θ)
    if m<0:
        Y = sqrt(2) * (-1)**m * Y.imag
    elif m>0:
        Y = sqrt(2) * (-1)**m * Y.real
    return Y

def R(r,n=1,l=0,Znuclear=1,a0=5.29177e-1):
    Zn = 2*Znuclear/n/a0
    ρ  = 2*Znuclear*r/(n*a0)
    N  = 1/factorial(2*l+1)*sqrt( Zn*Zn*Zn*factorial(n+l)/(2*n*factorial(n-l-1)) )
    R  = ρ**l*exp(-ρ/2)*N*hyp1f1(l+1-n,2*l+2,ρ)
    # https://en.wikipedia.org/wiki/Confluent_hypergeometric_function
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.hyp1f1.html
    return R

def Ψnlmxyz(n,l,m,X,Y,Z,X0=0,Y0=0,Z0=0):
    r = sqrt((X-X0)**2+(Y-Y0)**2+(Z-Z0)**2)
    θ = arccos((Z-Z0)/r)
    ϕ = arctan2((Y-Y0),(X-X0))+π
    return R(r,n=n,l=l)*realYlm(θ,ϕ,l=l,m=m)

def GetElementosMatriz(i,j,H2):
  ni,li,mi = 1,0,0
  sitioRi  = H2[i]
  Xi,Yi,Zi = sitioRi
  Ψnlm_i   = Ψnlmxyz(ni,li,mi,X,Y,Z,X0=Xi,Y0=Yi,Z0=Zi)
  nj,lj,mj = 1,0,0
  sitioRj  = H2[j]
  Xj,Yj,Zj = sitioRj
  Ψnlm_j   = Ψnlmxyz(nj,lj,mj,X,Y,Z,X0=Xj,Y0=Yj,Z0=Zj)

  DΨ_X  = gradient( Ψnlm_j,DX,axis=0,edge_order=2 )
  DDΨ_X = gradient(   DΨ_X,DX,axis=0,edge_order=2 )
  DΨ_Y  = gradient( Ψnlm_j,DY,axis=1,edge_order=2 )
  DDΨ_Y = gradient(   DΨ_Y,DY,axis=1,edge_order=2 )
  DΨ_Z  = gradient( Ψnlm_j,DZ,axis=2,edge_order=2 )
  DDΨ_Z = gradient(   DΨ_Z,DZ,axis=2,edge_order=2 )

  Kinetic = -ħ2_2m*conj(Ψnlm_i)*( DDΨ_X+DDΨ_Y+DDΨ_Z )
  Potencial = zeros_like(Kinetic)
  # 
  for n in [0,1]:
    Xp,Yp,Zp = H2[n]
    absDist  = sqrt( (X-Xp)*(X-Xp)+(Y-Yp)*(Y-Yp)++(Z-Zp)*(Z-Zp) )
    Potencial += e2_4πɛ0*conj(Ψnlm_i)*Ψnlm_j/absDist

  Hij = sum( (Kinetic+Potencial)*DX*DY*DZ )
  return Hij

def GetOverlap(i,j,H2):
  ni,li,mi = 1,0,0
  sitioRi  = H2[i]
  Xi,Yi,Zi = sitioRi
  Ψnlm_i   = Ψnlmxyz(ni,li,mi,X,Y,Z,X0=Xi,Y0=Yi,Z0=Zi)
  nj,lj,mj = 1,0,0
  sitioRj  = H2[j]
  Xj,Yj,Zj = sitioRj
  Ψnlm_j   = Ψnlmxyz(nj,lj,mj,X,Y,Z,X0=Xj,Y0=Yj,Z0=Zj)
  Sij = sum( (Ψnlm_i*Ψnlm_j)*DX*DY*DZ )
  return Sij



# Malla
N = 350*1J
xmin,xmax = -10,10
ymin,ymax = -10,10
zmin,zmax = -10,10
X,Y,Z = mgrid[xmin:xmax:N,ymin:ymax:N,zmin:zmax:N]
DX = X[(1,0,0)]-X[(0,0,0)]
DY = Y[(0,1,0)]-Y[(0,0,0)]
DZ = Z[(0,0,1)]-Z[(0,0,0)]


# Geometría
H2 = { 0:array([-0.50000, 0.00000,0.00000]),
       1:array([+0.50000, 0.00000,0.00000])}
# """

# Geometría
H2 = { 0:array([-0.50000, 0.00000,0.00000]),
       1:array([+0.50000, 0.00000,0.00000])}
h11 = GetElementosMatriz(0,0,H2)
h21 = GetElementosMatriz(1,0,H2)
h12 = GetElementosMatriz(0,1,H2)
h22 = GetElementosMatriz(1,1,H2)
print("h11 = {0}".format(h11.real))
print("h21 = {0}".format(h21.real))
print("h12 = {0}".format(h12.real))
print("h22 = {0}".format(h22.real))

S00 = GetOverlap(0,0,H2)
S10 = GetOverlap(1,0,H2)
S01 = GetOverlap(0,1,H2)
S11 = GetOverlap(1,1,H2)
print("S00 = {0}".format(S00.real))
print("S10 = {0}".format(S10.real))
print("S01 = {0}".format(S01.real))
print("S11 = {0}".format(S11.real))

h11 = -27.121729496271193
h21 = -20.255341653731666
h12 = -20.25534165373168
h22 = -27.12172949627123

S00 = 0.9999959401743611
S10 = 0.6165532580396124
S01 = 0.6165532580396124
S11 = 0.9999959401743653

2 Teorema de Bloch

Las soluciones de Equation 1 están dadas por el Teorema de Bloch.

Antes de enunciar el Teorema generemos una intuación sobre las soluciones de Equation 1. Grafiquemos quien es el término, recordando que tenemos una estructura cristalina

\[ \begin{align} V(\vec{r})=-\sum_{j=1}^{N} \frac{1}{4\pi\epsilon_0} \frac{eq_{\text{core},j}}{|\vec{R}_j-\vec{r}|} \end{align} \]

Code

from pylab import *
import plotly.graph_objects as go

def GenerarRed(n_1,n_2,a_1,a_2):
  # -- Desarrollo técnico --
  # Malla de enteros
  N1,N2 = meshgrid(n_1,n_2)
  # Cambiar a columnas
  N1    = N1.flatten()
  N2    = N2.flatten()
  N1y2 = column_stack((N1,N2))
  # Agrupar los vectores de red
  avec = [a_1,a_2]

  # Tomar todas las combinaciones
  Rx,Ry = dot( N1y2,avec ).T
  return Rx,Ry

a_1 = [1,0]
a_2 = [1/2,sqrt(3)/2]

n_1 = arange(-100,100)
n_2 = arange(-100,100)

Rx,Ry = GenerarRed(n_1,n_2,a_1,a_2)

sitio1x = Rx + (1/3)*( a_1[0]+a_2[0] )
sitio1y = Ry + (1/3)*( a_1[1]+a_2[1] )

sitio2x = Rx + (2/3)*( a_1[0]+a_2[0] )
sitio2y = Ry + (2/3)*( a_1[1]+a_2[1] )

npts = 60
xmin,xmax = -1,2
ymin,ymax = -1,2
X,Y = meshgrid(linspace(xmin,xmax,npts),linspace(xmin,xmax,npts))

V = zeros_like(X)
for n in range(len(Rx)):
    V += -1/(sqrt( (sitio1x[n]-X)*(sitio1x[n]-X) + (sitio1y[n]-Y)*(sitio1y[n]-Y) ))
    V += -1/(sqrt( (sitio2x[n]-X)*(sitio2x[n]-X) + (sitio2y[n]-Y)*(sitio2y[n]-Y) ))

Vmax = V.max()
fac  = 1.002

PLOTS = [go.Surface(x=X,y=Y,z=V,
                    cmin=fac*Vmax,
                    cmax=Vmax,
                    opacity=0.8,
                    showscale=False),
       go.Scatter3d(x=sitio1x,y=sitio1y,z=Vmax*ones_like(Rx),
                    mode='markers',showlegend=False),
       go.Scatter3d(x=sitio2x,y=sitio2y,z=Vmax*ones_like(Rx),
                    mode='markers',showlegend=False),
       go.Scatter3d(x=[0,a_1[0]],y=[0,a_1[1]],z=[Vmax,Vmax],
                    mode='lines',
                    line=dict(color='black',width=5),
                    showlegend=False),
       go.Scatter3d(x=[0,a_2[0]],y=[0,a_2[1]],z=[Vmax,Vmax],
                    mode='lines',
                    line=dict(color='black',width=5),
                    showlegend=False),
       go.Scatter3d(x=[a_1[0],a_1[0]+a_2[0]],y=[a_1[1],a_1[1]+a_2[1]],z=[Vmax,Vmax],
                    mode='lines',
                    line=dict(color='black',width=5),
                    showlegend=False),
       go.Scatter3d(x=[a_2[0],a_1[0]+a_2[0]],y=[a_2[1],a_1[1]+a_2[1]],z=[Vmax,Vmax],
                    mode='lines',
                    line=dict(color='black',width=5),
                    showlegend=False)]

fig = go.Figure( data=PLOTS )


camera = dict(
    eye=dict(x=0.5, y=2, z=0.1)
)

fig.update_layout(
    width=500, height=500,
    scene = dict(
        xaxis_title = 'x',
        yaxis_title = 'y',
        zaxis_title = 'V',
        xaxis = dict(range=[xmin,xmax],showticklabels=False),
        yaxis = dict(range=[ymin,ymax],showticklabels=False),
        zaxis = dict(nticks=2,range=[fac*Vmax,Vmax],showticklabels=False)
    ),
    scene_camera=camera
)

fig.show()

Dado que el potencial Equation 1 es periódico, se tiene el mismo entorno para todo punto, podemos esperar que la probabilidad de encontrar al electrón en un punto de la celda unitaria es la misma que el de una celda vecina.

Teorema de Bloch

Sea \(\hat{H}\) un Hamiltoniano de un electrón con un potencial periódico,

\[ \begin{align} \hat{H} = -\frac{\hbar^2}{2m}\hat{\nabla}^2 + U(\vec{r})\,,\quad \text{con}\;\; U(\vec{r})=U(\vec{r}+\vec{R}_{[\vec{n}]})\,, \end{align} \]

los eigenestados \(\psi\) de dicho Hamiltoniano pueden elegirse de la forma de una onda plana por una función con la periodicidad de la red. Por lo que podemos escribir las eigenfunciones de la siguiente forma,

\[ \begin{align} \psi_{n\vec{k}}(\vec{r}) = e^{i \vec{k} \cdot\vec{r} } u_{n\vec{k}}(\vec{r})\,, \end{align} \]

donde

\[ \begin{align} u_{n\vec{k}}(\vec{r}+\vec{R}_{[\vec{n}]}) = u_{n\vec{k}}(\vec{r})\,,\quad \forall \,\vec{R}_{[\vec{n}]}\,. \end{align} \]

\[ z = |z|e^{i\theta} \]

Notamos que:

1. \(\psi_{\vec{k}}(\vec{r})\) no es periódica excepto si \(\vec{k}=0\)

2. Si \(\vec{k}\notin {\rm 1BZ}\) entonces \(\vec{k}=\vec{G}+\vec{k}_{\text{1BZ}}\) donde \(\vec{G}\) es un vector de la red recíproca.

3. \[ \begin{align} \psi_{\vec{k}}(\vec{r}+\vec{R}) &= u_{\vec{k}}(\vec{r}+\vec{R}) e^{i\vec{k}\cdot (\vec{r}+\vec{R})} \\ &= u_{\vec{k}}(\vec{r}) e^{i\vec{k}\cdot \vec{r}} e^{i\vec{k}\cdot \vec{R}}\\ &=\psi(\vec{r}) e^{i\vec{k}\cdot \vec{R}} \end{align} \]

Alternativamente, el Teorema de Bloch nos dice que la solución adquiere una fase al trasladarse en el espacio.

2.1 Propiedades ópticas

¿Cuál es la longitud de onda de la luz absorbida, transmitida o reflejada?

Luz visible 740 nm (1.7 eV) red 400 nm violet (3.1 eV)

E gap= min conduction band energy max valance band energy

Filled valance band and empty conduction band

Aislante con Egap > 3.1 eV es transparente Diamond, quartz, glass, alumina, gallium nitride

Samller band gap, only photons con \(\hbar \omega>E_{\rm GAP}\) pueden absorberse

Colores

GAP indirecto y GAP indirecto

Direc GAP transitions -> strong light absortion

Indirect transitions -> weak absortion

1 photon -> 1 electron (excited) + 1 photon (excited)

Imperfect momentum conservation From Disorder / impurities

“Cartoon”

Cartoon

Code

import numpy as np
import matplotlib.pyplot as plt

# wavelength range, nm
lmin, lmax = 250, 1000
x = np.linspace(lmin, lmax, 1000)
# A wave with a smoothly increasing wavelength
wv = (np.sin(10 * np.pi * x / (lmax+lmin-x)))[::-1]

fig = plt.figure()
ax = fig.add_subplot(111, facecolor='k')
ax.plot(x, wv, c='w', lw=2)
ax.set_xlim(250,1000)
ax.set_ylim(-2,2)

# Label and delimit the different regions of the electromagnetic spectrum
ax.text(310, 1.5, 'UV', color='w', fontdict={'fontsize': 20})
ax.text(530, 1.5, 'Visible', color='k', fontdict={'fontsize': 20})
ax.annotate('', (400, 1.3), (750, 1.3), arrowprops={'arrowstyle': '<|-|>',
                                                    'color': 'w', 'lw': 2})
ax.text(860, 1.5, 'IR', color='w', fontdict={'fontsize': 20})
ax.axvline(400, -2, 2, c='w', ls='--')
ax.axvline(750, -2, 2, c='w', ls='--')
# Horizontal "axis" across the centre of the wave
ax.axhline(c='w')
# Ditch the y-axis ticks and labels; label the x-axis
ax.yaxis.set_visible(False)
ax.set_xlabel(r'$\lambda\;/\mathrm{nm}$')

# Finally, add some colourful rectangles representing a rainbow in the
# visible region of the spectrum.
# Dictionary mapping of wavelength regions (nm) to approximate RGB values
rainbow_rgb = { (400, 440): '#8b00ff', (440, 460): '#4b0082',
                (460, 500): '#0000ff', (500, 570): '#00ff00',
                (570, 590): '#ffff00', (590, 620): '#ff7f00',
                (620, 750): '#ff0000'}
for wv_range, rgb in rainbow_rgb.items():
    ax.axvspan(*wv_range, color=rgb, ec='none', alpha=1)
plt.show()

Silicio

Details matter

1. Surface
2. Impurities

Estados entre los GAPS

Diamond

Dopant = foreing chemical otherwise pure substance to change its properties

Mind the GAP